setwd('/home/ruser/zychw')

current_libs <- .libPaths()
your_personal_lib_path <- "~/R/library"
.libPaths(c(current_libs, your_personal_lib_path))


# 下载 TCGA-GBM STAR-Counts数据 -----------------------------------------------

library('TCGAbiolinks')
## 设置下载信息
CancerProject <- "TCGA-GBM"
query <- GDCquery(project = CancerProject, 
                  data.category = "Transcriptome Profiling", 
                  data.type = "Gene Expression Quantification", 
                  workflow.type = "STAR - Counts")
## 开始下载
GDCdownload(query) 
dataGBM <- GDCprepare(query)

dataGBM <- dataGBM[,!is.na(dataGBM$paper_IDH.status)]


# 使用 DESeq2 进行差异表达分析 ------------------------------------------------------

library('DESeq2')
## 1.得到paper_IDH.status分组
ddsSE <- DESeqDataSet(dataGBM, design = ~ paper_IDH.status)
## 2.过滤低表达数据
keep <- rowSums(counts(ddsSE)) >= 10 
ddsSE <- ddsSE[keep,]
## 3.使用DESeq2进行差异表达分析
ddsSE <- DESeq(ddsSE)

## 展示结果
# results(ddsSE)
# 原始counts数据在 ddsSE@assays@data$counts



# 添加注释信息 ------------------------------------------------------------------

library('org.Hs.eg.db')
library('stringr') 

## 将 ENSEMBL 后的小数点删除
gene_ids <- str_split(rownames(ddsSE),'[.]',simplify = T)[,1]
rownames(ddsSE) <- gene_ids

## 查询得需要注释信息表格
annot <- select(org.Hs.eg.db, keys = gene_ids,
                column = c('ENSEMBL','ENTREZID','GENENAME','GENETYPE','SYMBOL'),
                keytype = 'ENSEMBL')

# 记录数据至 deseq2_results.txt ----------------------------------------------------

## 将差异表达数据写入 deseq2_result.txt
## 取出 IDH WT vs mutant 的数据
de <- results(object = ddsSE, contrast=c("paper_IDH.status", "WT", "Mutant"))
## 按 ENSEMBL==ID 连接数据表格与注释表格
de_symbols <- merge(data.frame(ID=rownames(de), de, check.names=FALSE), annot, 
                    by.x="ID", by.y="ENSEMBL", all=F)
## 文件记录
write.table(de_symbols, "deseq2_results.txt", quote=F, col.names=T, row.names=F, sep="\t")


# 筛选差异表达的基因 ---------------------------------------------------------------

# ## 1.log2FoldChange

# ### 计算 log2FoldChange
# norm_counts <- log2(counts(ddsSE, normalized = TRUE)+1)
# ### 连表注释
# norm_counts_symbols <- merge(data.frame(ID=rownames(norm_counts), norm_counts, check.names=FALSE), 
#                              annot, by.x="ID", by.y="ENSEMBL", all=F)
# ### 文件记录 
# write.table(norm_counts_symbols, "normalized_counts_log2_star.txt", quote=F, col.names=T, row.names=F, sep="\t")
# 
# ## 2.padj < 0.05

# ### 筛选 padj < 0.05 的基因
# de_select <- de_symbols[de_symbols$padj < 0.05 & !is.na(de_symbols$padj),]
# ### 文件记录 
# write.table(de_select, "deseq2_selection_padj005.txt", quote=F, col.names=T, row.names=F, sep="\t")


# 画火山图 --------------------------------------------------------------------

library('magrittr')
library('ggplot2')

# deg_deseq2_ret <- readr::read_delim("deseq2_results.txt")
deg_deseq2_ret <- de_symbols

## 将所有标记为 "NOTSIG"
deg_deseq2_ret$diffexpressed <- "NOTSIG"
## 再将 log2Fold > 1 且 padj < 0.05 标记为 "UP" 
deg_deseq2_ret$diffexpressed[deg_deseq2_ret$log2FoldChange > 1 & 
                               deg_deseq2_ret$padj < .05] <- "UP"
## 将 log2Fold < -1 且 padj < 0.05 标记为 "DOWM" 
deg_deseq2_ret$diffexpressed[deg_deseq2_ret$log2FoldChange < -1 & 
                               deg_deseq2_ret$padj < .05] <- "DOWN"
## 设置颜色
mycolors <- setNames(c("cornflowerblue","grey", "firebrick"), c("DOWN","NOTSIG", "UP"))


library('ggrepel')
## 按log2FoldChange降序排序，取前10
top_10_rows <- deg_deseq2_ret[order(abs(deg_deseq2_ret$log2FoldChange), decreasing = TRUE),][1:10,]


## 绘图
deg_deseq2_ret %>% 
  ggplot( aes(x=log2FoldChange, y=-log10(padj), 
              color=diffexpressed)) + 
  geom_point() +    
  geom_text_repel(data = top_10_rows, aes(label=SYMBOL)) + 
  geom_vline(xintercept=c(-1,1), col="red", linetype=2) +
  geom_hline(yintercept=-log10(.05), col="red", linetype=2) +
  theme_minimal(base_size = 16) +
  scale_color_manual(values = mycolors) -> p4

## 写入volcano.pdf文件
pdf("volcano.pdf")
p4
dev.off()


# 提取David所需的数据 ---------------------------------------------------------

## 得到 Entrez IDs
deg_deseq2_ret %>% 
  dplyr::filter(diffexpressed != "NOTSIG") %>% 
  dplyr::pull(ENTREZID) %>% 
  as.character() -> degs


degs %>% as.data.frame() %>% 
  readr::write_csv("toDavid.csv", col_names = F)
